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Abstract 

By using a powerful reductive perturbation technique, or a multiscale analysis, a generic 
Kadomtsev-Petviashvili evolution equation governing the propagation of femtosecond spatiotempo- 
ral optical solitons in quadratic nonlinear media beyond the slowly- varying envelope approximation 
is put forward. Direct numerical simulations show the formation, from adequately chosen few-cycle 
input pulses, of both stable line solitons (in the case of a quadratic medium with normal dispersion) 
and of stable lumps (for a quadratic medium with anomalous dispersion). Besides, a typical exam- 
ple of the decay of the perturbed unstable line soliton into stable lumps for a quadratic nonlinear 
medium with anomalous dispersion is also given. 

PACS numbers: 42.65.Tg, 42.65.Rc, 05.45.Yv 
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I. INTRODUCTION 



Over the past two decades, the development of colliding-pulse mode-locked lasers and self- 
started solid-state lasers has been crowned with the generation of few-cycle pulses (FCPs) 
'or wavelengths from the visible to the near infrared (for a comprehensive review see Ref. 
l|). Presently, there is still a great deal of interest in the study of propagation character- 
istics of these FCPs in both linear and nonlinear media. Thus the comprehensive study of 
intense optical pulses has opened the door to a series of applications in various fields such as 
light matter interaction^ high-order harmonic generation, extreme nonlinear optics [2|, and 



sics Q,0]. 



attosecond physics [3|, |4|. The theoretical studies of the physics of FCPs concentrated on 
three classes of main dynamical models: (i) the quantum approach 5|-|8| , (ii) the refinements 
within the framework of the slowly varying envelope approximation (SVEA) of the nonlin- 
ear Schrodinger-type envelope equations |9h12|. and non-SVEA models 13Ml7j|. In media 
with cubic optical nonlinearity (Kerr media) the physics of (l+l)-dimensional FCPs can be 



adequately described beyond the SV 
modified Korteweg-de Vries (mKdV) 







n 



A by using different dynamical models, such as the 



13j, sine-Gordon (sG) [1^, ll5(] , or mKdV-sG equations 



17| . The above mentioned partial differential evolution equations admit breather solu- 
tions, which are suitable for describing the physics of few-optical-cycle solitons. Moreover, 



a non-integrable generalized Kadomtsev-Petviashvili (KP) equation [18[ (a two-dimensional 
version of the mKdV model) was also put forward for describing the (2+1) -dimensional 
few-optical-cycle spatiotemporal soliton propagation in cubic nonlinear media beyond the 



SVEA 



19 



20]. 



Notice that of particular interest for the physics of the (l+l)-dimensional FCPs is the 
mKdV-sG equation; thus by using a system of two-level atoms, it has been shown in Ref. 



161 ] that the propagation of ultrashort pulses in Kerr optical media is fairly well described by 



a generic mKdV-sG equation. It is worthy to mention that this quite general model was also 
derived and studied in Refs. 



21 



221 ] . However, another model equation for describing the 



physics of FCPs, known as the short-pulse equation (SPE) has been introduced, too [231-126] 



and some vectorial versions of SPE have been also investigated 



27 



29j. Other kind of SPE 



containing additional dispersion term has 
version of the SPE was also put forward 
been studied, too 



aeen introduced in Ref. I30J1. A multi-dimensional 



311 ] and self-focusing and pulse compression have 



32| . Notice that the SPE model has been considered again in its vectorial 
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version and the pulse self- compression and FCP soliton propagation have been investigated 



in detail 



33] . The main dynamical (l+l)-dimensional FCP models mentioned above have 



been revisited recently 34J; it was thus proved that the generic dynamical model based on 
the mKdV-sG partial differential equation was able to retrieve the results reported so far in 
the literature, and so demonstrating its remarkable mathematical capabilities in describing 
the physics of (l+l)-dimensional FCP optical solitons 34]. 

As concerning the extension of these studies to other relevant physical settings, it was 
shown recently that a FCP launched in a quadratically nonlinear medium may result in the 
formation of a (l+l)-dimensional half-cycle soliton (with a single hump) and without any 
oscilatting tails 35] . It was proved that the FCP soliton propagation in a quadratic medium 
can be adequately described by a KdV equation and not by mKdV equation. Notice that in 
Ref. 35[ it was considered a quadratic nonlinearity for a single wave (frequency), and that no 



effective third order nonlinearity was involved (due to cascaded second order nonlinearities). 
This is in sharp contrast with the nonlinear propagation of standard quadratic solitons (alias 
two-color solitons) within the SVEA where two different waves (frequencies) are involved, 
namely a fundamental frequency and a second harmonic 36|- 43]. 

As mentioned above, there are only a few works devoted to the study of (2+l)-dimensional 
few-optical-cycle solitons, where an additional spatial transverse dimension is incorporated 



into the model; see e.g., Refs. 



15 



19 



20(, where the propagation of few-cycle optical 



pulses in a collection of two-level atoms was investigated beyond the traditional SVEA. The 
numerical simulations have shown that in certain conditions, a femtosecond pulse can evolve 



into a stable few-optical cycle spatiotemporal soliton 15 



19 



The aim of this paper is to derive a generic partial differential equation describing the 
dynamics of (2+l)-dimensional spatiotemporal solitons in quadratic nonlinear media beyond 
the SVEA model equations. Our study rely on the previous work 35[ where the theory of 
half-cycle (1+1) -dimensional optical solitons in quadratic nonlinear media was put forward. 
This paper is organized as follows. In Sec. II we derive the generic KP equation governing 
the propagation of femtosecond spatiotemporal solitons in quadratic nonlinear media beyond 
the SVEA. To this aim we use the powerful reductive perturbation technique up to sixth 
order in a small parameter e. We have arrived at the conclusion that when the resonance 
frequency is well above the inverse of the typical pulse width, which is of the order of a 
few femtoseconds, the long-wave approximation leads to generic KP I and KP II evolution 
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equations. These evolution equations and their known analytical solutions, such as line 
solitons and lumps are briefly discussed in Sec. III. Direct numerical simulations of the 
governing equations are given in Sec. IV, where it is shown the generation of both stable 
line solitons (for the KP II equation) and of stable lumps (for the KP I equation), and a 
typical example of the decay of the perturbed unstable line soliton of the KP I equation into 
stable lumps is also provided. Section V presents our conclusions. 



II. DERIVATION OF THE KP EQUATION 



We consider the same classical and quantum models developed in Ref. 35[, but with an 
additional transverse spatial dimension (we call it the y-axis). In the quantum mechani- 
cal approach, the medium composed by two-level atoms is treated using the density-matrix 
formalism and the dynamics of electromagnetic field is governed by the Maxwell- Bio ch equa- 
tions. Thus the evolution of the electric field component E along the polarization direction 
of the wave (the x-axis) is described by the (2+l)-dimensional Maxwell equations, which 
reduce to 

(St + d 2 z ) E = i# (E + 4vrP) , (1) 

where P is the polarization density. 

a. Classical model. First we consider a classical model of an elastically bounded elec- 
tron oscillating along the polarization direction of the wave. The wave itself propagates in 
the ^-direction. The evolution of the position x of an atomic electron is described by an 
anharmonic oscillator (here damping is neglected) 

d -^ + tt 2 x + ax 2 = —E. (2) 
at 1 m 

The parameter a measures the strength of the quadratic nonlinearity. The polarization 
density is P = —Nex, where N is the density of atoms. Here Q is the resonance frequency, 
(— e) is the electron charge, m is its mass, and E is the electric field component along the 
x-axis. 

b. Quantum model. We consider a set of two- level atoms with the Hamiltonian 



(3) 
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where f2 = ojb — uj a > is the frequency of the transition. The evolution of the electric 
field E is described by the wave equation ([T]). The light propagation is coupled with the 
medium by means of a dipolar electric momentum \i directed along the same direction x as 
the electric field, according to 

H = H - pE, (4) 
and the polarization density P along the ^-direction is 

P = NTr(pp), (5) 

where N is the volume density of atoms, and p is the density matrix. The density-matrix 
evolution equation (Schrodinger equation) is written as 

ihd tP =[H,p]+K, (6) 

where 1Z is a phenomenological relaxation term. In the case of cubic nonlinearity, it has 
been shown that it was negligible [ljj], and it will be neglected here as well. 
In the absence of permanent dipolar momentum, 



f=\ (7) 




is off-diagonal, and the quadratic nonlinearity is zero. 

The quadratic nonlinearity can be phenomenologically accounted for as follows: it corre- 
sponds to a deformation of the electronic cloud induced by the field E, hence to a dependency 
of the energy of the excited level b with respect to E, i.e., to a Stark effect. The nonlinearity 
is quadratic if this dependency is linear, as follows 

(j b — > Ub — aE. (8) 

This contribution can be included phenomenologically in the Maxwell-Bloch equations by 
replacing the free Hamiltonian H with 



(9) 
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c. Scaling. Transparency implies that the characteristic frequency u w of the considered 
radiation (in the optical range) strongly differs from the resonance frequency Q of the atoms, 
hence it can be much higher or much lower. We consider here the latter case, i.e., we assume 
that oj w is much smaller than Q. This motivates the introduction of the slow variables 



e it 



z 
V 



e z. 



n 



(10) 



e being a small parameter. The delayed time r involves propagation at some speed V to be 
determined. It is assumed to vary slowly in time, according to the assumption u w <C Q. The 
pulse shape described by the variable r is expected to evolve slowly in time, the corresponding 
scale being that of variable (. The transverse spatial variable y has an intermediate scale as 
usual in KP-type expansions 44]. 

A weak amplitude assumption is needed in order that the nonlinear effects arise at the 
same propagation distance scale as the dispersion does. Next we use the reductive perturba- 



tion method as developed in Ref. 
series of a small parameter e: 



44J . To this aim we expand the electric field E as power 



E = E 2 E 2 4~ s^E^ 



e A E 4 + ..., 



(11) 



as in the standard KdV-type expansions 44J. The polarization density P is expanded in the 
same way. 



d. Order by order resolution The computation is exactly the same as in 35| up to order 
e 4 in the Schrodinger equation (|6]) or classical dynamical equation fl2]), and up to order e 5 
in the wave equation (CQ). The fourth order polarization density P4 has the same expression 



as m 



35j | . Then the wave equation at order e 6 yields the evolution equation 

d c d T E 2 = AdtE 2 + Bd 2 T (E 2 ) 2 + ^d 2 v E 2 , 

which is a generic KP equation. 

The coefficients A and B in the above (2+l)-dimensional evolution equation are 



(12) 



A 



AnVNe 2 



2mc 2 n 4 

in the classical model, and respectively, 

4vriV|/i| 2 



B 



A 



nchQ 3 



B 



-4nVaNe 3 
2m 2 c 2 Q 6 



-AirNa\n\ 
nch 2 Q 2 



(13) 



(14) 



in the corresponding quantum model. 

Notice that the dispersion coefficient A and the nonlinear coefficient B can be written in 
a general form as 



.4 



1 d 3 k 



B 



6 du 3 
-2tt 



1 d 2 n 



nc 



2c du 2 



oj=0 
(2), 



cv=0 



uj=0 



(15) 



(16) 



where n is the refractive index of the medium and u is the wave pulsation (for a detailed 
discussion of the corresponding expressions of the dispersion and nonlinear coefficients in 
the case of cubic nonlinear media see, e.g., Ref. Ijj]). The velocity is obviously V = c/n. 



III. TWO KP EQUATIONS: KP I AND KP II 

The KP equation ( 1121) is reduced to 

dzdru = d^u + d^u 2 + adyU, (17) 
where a = sgn(A), by the change of variables 

(V A 1 / 3 

C = Z, r = A 1 / 3 T, V =\A\ 1 ^-Y, E 2 = —u. (18) 

For our starting models (both the classical and quantum ones), the dispersion coefficient 
A is positive. It is the case of the normal dispersion. Hence a — +1 and Eq. f|T7|) is the 
so-called KP II equation. Both KP I and KP II are completely integrable by means of the 



inverse scattering transform (1ST) method 45], however the mathematical properties of the 



solutions differ. KP II admits stable line solitons, but not stable localized soliton solutions. 
The line soliton has the analytical expression 



k,2 = — sech ^ p 



r + ar]^ I ^ + 4p 2 A ) ( 



(19) 



where a and p are arbitrary constants and A, B have the expressions given in the previous 
Section. 

If Eq. (|T2|) is generalized to other situations using the expressions ( |T5|) -( TT6|) for the 
dispersion and nonlinear coefficients, and if we assume an anomalous dispersion, i.e., A < 0, 
then a = — 1 and Eq. ( ITTl) becomes the so-called KP I equation, for which the above line 
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FIG. 1: (Color online) The lump solution of the KP I equation. Here the parameters are A = — 1, 
B = — 1, V = 1, p r = — 1, and pi = 2. 



soliton is not stable, but which admits stable lumps [46|, |47] as 

„ _ 12 A [- (t' + p r yf+p 2 y' 2 + 3/p 2 ] 
B [(t'+p r y') 2 +p 2 y' 2 + 3/p 2 ] 

with 

t' = T+{p 2 r +p 2 )A(, (21) 



-2A 

y' = V\l—-2PrA(, (22) 

where p r and Pi are arbitrary real constants. The analytical lump solution written above is 
plotted in Fig. [TJ 

The size of the lump can be evaluated as follows. Assume for the sake of simplicity that 
the lump exactly propagates along z, i.e., that p r = 0. Then Eq. ( 1201) has the form 

/ J./2 ,2 \ - 2 

E 2 « ( TTT + TT-Z + 1 > ( 23 ) 



which gives the normalized transverse width and duration as w;^ = V^/p 2 and = v3/Pi> 
respectively, while the scaled peak amplitude of the lump is E-imax = ^Ap 2 /B. Using Eqs. 
(12~TT)-(122"1) and the scaling ( |T0l) . we get the corresponding width and duration in physical 
units, as 

yg / V vg 



^ = :xW^77' U7 * = — ' ( 24 ) 



S 



while the peak amplitude is E v 
cients A and B, we obtain 



£ 2 E2max • Using expressions (jl~5l) and (|T6|) of the coeffi- 



w 2 c 



-3nn'' 



(25) 



\J-3nn" ' m 7rx (2) Wt' 

where n" and x are shortcuts for d 2 n/ duo 2 \ uJ=0 and x < - 2 - ) (2o;; w, w) | =0 , respectively. Hence 
the transverse width of the lump decreases as the pulse duration does and the peak amplitude 
increases. However, the initial assumptions and scaling ( flOj) imply that the field varies much 
slower in the y-direction than longitudinally, hence the width of the lump should, in principle, 
be large with respect to its length. If this assumption is not satisfied, the validity of the KP 
model is questionable. 



IV. NUMERICAL SIMULATIONS OF KP I AND KP II EQUATIONS 



The KP equation (ITT)) is solved by means of the fourth order Runge Kutta exponential 
time differencing (RK4ETD) scheme 48| . It involves one integration with respect to r. 
The inverse derivative is computed by means of a Fourier transform, which implies that the 
integration constant is fixed so that the mean value of the inverse derivative is zero, but also 
that the linear term is replaced with zero, i.e., the mean value of the function d 2 u is set to 
zero. For low frequencies, the coefficients of the RK4ETD scheme are computed by means 
of series expansions, to avoid catastrophic consequences of limited numerical accuracy. 



A. KP II 



For a normal dispersion (n" > 0), as in the case of the two-level model above, Eq. ([12 
is KP II. Then, starting from an input in the form of a Gaussian plane wave packet 

\2 



u = u exp 



T-Tx 



wr 



cos 



y{r-7i: 



where 



^ -A / Y 2 
T\ = —— exp 



IV; 



(26) 



(27) 



accounts for a transverse initial perturbation, we have checked numerically (see Fig. 2) 
that (i) The FCP input transforms into a half-cycle pulse, line soliton solution of KP II 
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equation, whose temporal profile is the sech-square shaped KdV soliton, plus the accompa- 
nying dispersing- diffracting waves, and that (ii) The initial perturbation vanishes and the 
line soliton becomes straight again during propagation. 



According to the stability analysis 45|, the stabilization of the perturbed line-soliton 
occurs if the perturbation is not too large. Hence the essential requirement to evidentiate 
the above features is the one- dimensional problem of the formation of a KdV soliton from 
a Gaussian input. In principle, this problem is solved by the 1ST: a wide pulse with high 
enough intensity is expected to evolve into a finite number of solitons, plus some dispersive 
waves, called 'radiation' in the frame of the 1ST method. This has been addressed numeri- 



cally in [35|: it has been seen that, for a FCP-type input, the amount of radiation was quite 
important, and the exact number and size of solitons might depend on the carrier-envelope 
phase. 

Numerically, the finite size of the computation box and periodic boundary conditions 
prevent the dispersive wave from being totally spread out. Hence, in order to evidence 
numerically the phenomenon mentioned above, we need to reduce the energy of the dispersive 
part to the lowest level possible, and produce a single soliton. Very roughly, one can say 
that each positive half-cycle oscillation with high enough amplitude results in one soliton. 
The initial carrier envelope phase should be close to zero, the main oscillation have a large 
enough amplitude, and the pulse be short enough so that only one peak is emitted. On the 
other hand, the length of the KdV soliton (fT9l) is 2/p, while a half-cycle of the pulse ( 1261) 
obviously has duration A/2. On the other hand, the amplitude of the KdV soliton ( |T9l) is 
Qp 2 A/B, to be compared to the amplitude uq of pulse (1261) . In a first approximation, the 
half-cycle matches the KdV soliton if both the durations and amplitude are the same, hence 
if 

"° = Pb- (28) 

Numerical computation confirms that relation (|28|) between the wavelength and amplitude 
of the Gaussian pulse is roughly optimal. 

Recall that our numerical computations are performed using the normalized form ( II 7p of 
the KP equation, which corresponds to A — B — 1. We observed that, for the particular 
values of A and Wt pertaining to Fig. [2, for no > 5, the soliton emerges from the dispersive 
wave, while, at optimum ( |28|) (m — 11), the soliton amplitude is about three times the 
amplitude of the dispersive wave. 
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FIG. 2: (Color online) Shape and profile of the line soliton. a-b: input data, c-d: after propagation 
Z = 15 (exactly 14.95). Parameters are A = 3, uq = 10.7, wt = 1.4, and w y = 3.87. 



A stable plane-wave half-cycle quadratic optical soliton can thus be formed from a FCP 
input. Further, the quadratic nonlinearity may restore the spatial coherence of the soliton. 

B. KP I 

Let us now assume a quadratic nonlinear medium with anomalous dispersion (n" < 0), in 
which the KP equation ffl2|) is valid. It reduce to Eq. (fT7|) with a — — 1, i.e., KP I equation. 
We consider an initial data of the form 

n = n ex P (-^-^)cos(|T), (29) 

with arbitrary values of the initial amplitude uq, width w y , and duration wt of the pulse. The 
evolution of the input FCP into a lump is shown on Fig. 3; notice the remaining dispersive 
waves. The fact that they seem to be confined in the ^-direction is a numerical artifact, due 
to the fact that absorbing boundary conditions have been added in the ^-direction only. 
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FIG. 3: (Color online) Shape and profile of the lump, a-b: input data, c-d: after propagation 
Z = 15 (exactly 14.95). Parameters are A = 3, uq = 6.12, wt = 1.2, and w y = 6.36. 




(a) (b) 

FIG. 4: (Color online) Generation of lumps from a perturbed unstable line soliton. a : input data, 
b: after propagation Z = 1.45. Parameters are A = 3, uq = 10.7, wt = 1.4, and L = 14.14, 



Hence the diffracting part of the dispersive waves rapidly decreases as if the width of the 
numerical box were infinite, while the dispersive waves which propagate along the Z-axis 
remain confined in the numerical box (with length 160), and hence decays slower as it would 
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happen for a single pulse. 

Figure 4 shows the decay of a perturbed unstable line soliton of KP I equation into lumps. 
The input is given by Eq. (I26p . with the transverse perturbation 



where the wavelength L of the transverse perturbation is chosen so that 2L = 28.3 is the 
width of the numerical box, in accordance with the periodic boundary conditions used in 
our numerical computation. This transversely perturbed line-soliton, in contrast with the 
case of KP II for normal dispersion shown in Fig. 2, does not recover its initial straight 
line shape and transverse coherence, but breaks up into localized lumps. Notice that the 
conditions on the parameters needed to obtain a result comparable to that shown in Fig. 4 
are much less restrictive than the conditions required for the numerical observation of line 
solitons; see Fig. 2 and the detailed discussion given above on the optimal conditions of 
generation of line solitons of KP II equation. Hence spatiotemporal quadratic solitons may 
form from a FCP input, their transverse focusing follows from any transverse perturbation 
of the incident plane wave. 

V. CONCLUSION 

In conclusion, we have introduced a model beyond the slowly varying envelope approxi- 
mation of the nonlinear Schrodinger-type evolution equations, for describing the propagation 
of (2+l)-dimensional spatiotemporal ultrashort optical solitons in quadratic nonlinear me- 
dia. Our approach is based on the Maxwell-Bloch equations for an ensemble of two level 
atoms. We have used the multiscale approach up to the sixth-order in a certain small per- 
turbation parameter and as a result of this powerful reductive perturbation method, two 
generic partial differential evolutions equations were put forward, namely, (i) the Kadomtsev- 
Petviashvili I equation (for anomalous dispersion), and (ii) the Kadomtsev-Petviashvili II 
equation (for normal dispersion). Direct numerical simulations of these governing partial 
differential equations show the generation of both stable line solitons (for the Kadomtsev- 
Petviashvili II equation) and of stable lumps (for the Kadomtsev-Petviashvili I equation). 
Moreover, a typical example of the decay of the perturbed unstable line solitons of the 
Kadomtsev-Petviashvili I equation into stable lumps is also given. 




(30) 
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The present sudy is restricted to (2+1) dimensions, however, due to the known properties 
of the Kadomtsev-Petviashvili equation in (3+1) dimensions |49|, |50|, analogous behavior is 
expected in three dimensions, i.e., stability of the spatial coherence of a plane wave few- 
cycle quadratic soliton for normal dispersion, and spontaneous formation of spatio-temporal 
few-cycle 'light bullets' in the case of anomalous dispersion. 
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